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ABSTRACT 

The  delineation  of  drainage  gullies,  in  mapping  operations, 
is  a  highly  labor  intensive  effort.  ETL  has  initiated  a 
preliminary  study  to  automate  this  manual  task.  This  has  led 
to  a  computer  algorithm  for  extracting  the  drainage  channels 
from  digital  elevation  data.  The  algorithm,  which  is  based 
on  the  theory  of  critical  numbers,  provides  both  the  location 
and  orientation  of  the  gullies.  This  paper  discusses  the 
algorithm  and  the  preliminary  results  from  three  test  cases. 


INTRODUCTION 

Rivers  and  streams  carry  the  water  from  some  68.7  percent  of 
the  earth's  landarea;  therefore,  the  total  number  of  drainage 
gullies  is  immense  (Gregory,  1973).  For  the  United  States 
alone,  there  are  over  five  million  kilometers  of  river 
channels.  This  figure  takes  into  consideration  only  those 
channels  which  would  be  considered  as  first  order  and  higher 
on  maps  of  a  scale  of  1/62,300  (Leet,  1971).  Currently,  the 
delineation  of  these  channels,  in  mapping  operat ions  ,  is  a 
highly  labor  intensive  task  and  automation  is  necessary  to 
decrease  the  production  time  and  cost. 

ETL  has  initiated  a  preliminary  study,  having  as  its  objective, 
the  sutomation  of  drainage  delineation  in  the  mapping  procesa. 
The  approach  uses  terrain  elevation  data,  such  as  those  which 
are  currently  produced  automatically  by  the  current  map 
production  systema.  This  elevation  data  is  processed  by  an 
algorithm,  which  extracts  both  location  and  direction  of 
drainage  gullies. 
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DELINEATION  METHOD 


The  computer  algorithm  for  the  delineation  of  drainage 
channnels  solely  from  digital  terrain  elevation  data  is  baaed 
on  an  idea  similiar  to  the  theory  of  critical  numbers.  In 
mathematics,  a  critical  number,  C,  of  a  continuous  function 
F  is  a  value  such  that  F(C)  is  a  relative  extremum  for  some 
closed  interval  [H,N]  and  either  H<C<B  or  M>C>  11  (Johnson, 
1970).  An  elevation  value  would  be  a  relative  extremum,  if 
for  some  given  neighborhood,  it  is  either  a  maximum  or  minimum 
value.  For  gully  delineation,  the  value  must  be  a  minimum; 
while,  for  ridge  delineation  only  the  maximums  are  considered. 

Given  digital  elevation  data  in  a  matrix  format,  then  by 

definition,  the  general  neighbornood  of  the  elevation  value 
is  its  eight  surrounding  values.  These  eight  values  form 
four  individual  neighborhoods  of  the  central  e levat ion 'value . 
In  Figure  1,  consider  E  as  the  elevation  value  and  points  1 
through  8  as  its  surrounding  values.  The  four  neighborhoods 
of  E  to  be  examined  are:  [1,8],  [2,7],  [3,6],  and  [4,5]. 

1  2  3 

4  E  5 
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Figure  1 


Each  of  the  neighborhoods  must  be  considered  since  a  drainage 
gully  can  have  any  one  or  a  combination  of  the  four  directional 
orientations.  If  the  test  results  prove  E  to  be  a  minimum 
for  any  of  the  cases,  then  the  appropriate  tag  is  recorded. 
This  tag  represents  both  the  location  and  orientation  of  the 
gully.  The  tagging  code,  which  has  been  implemented  in  this 
program,  is  a  binary  system.  A  "1"  is  noted  if  the  value  is 
a  minimum  for  the  neighborhood,  and  a  "0"  otherwise.  This 
scheme  was  organised  in  this  manner  to  eimplify  the  graphic 
subroutines.  Using  the  above  nine  element  matrix  as  an  example, 
the  four  rule'h  for  the  tagging  scheme  are: 


1  -  If  E  is  a  minimum  when 

then  the  tag  is  0001. 

2  -  If  E  is  a  minimum  when 

then  the  tag  is  0010. 

3  -  If  E  is  a  minimum  when 

then  the  tag  is  0100. 

4  -  If  E  is  a  minimum  when 

then  the  tag  is  1000. 
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Using  the  above  rules  on  the  element  ”422"  of  Figure  2,  the 
tag  would  be  0010.  This  results  since  it  was  in  the  comparison 
of  the  neighborhood  [425,425]  for  which  "422”  was  found  to 
be  a  minimum. 
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This  tagging  scheme  allows  for  the  previously  mentioned 
multiple-orientation  problem.  For  example,  given  the  following 
elevation  matrix  of  a  gully  with  a  two  directional  orientation, 


the  tagging  scheme,  step-by- 

step  , 

would  yield 
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Figure  3 

1  -  After  checking  the  neighborhood  (101,120)  the  tag  becomes 

0001  . 

2  -  After  checking  the  neighborhood  (83,85)  the  tag  remains 

0001  . 

3  -  After  checking  the  neighborhood  (85,130)  the  tag  remains 

0001  . 

4  -  After  checking  the  neighborhood  (120,100)  the  tag  becomes 

1001. 

Thus  the  resulting  tag  is  1001  -  a  multi-directional  orientation. 

In  order  to  process  large  matrices  of  elevation  data,  the 
program  first  retrieves  three  consecutive  columns  of  data. 
The  algorithm  then  does  the  testing  sequentially  on  3  by  3 
subset  matrices,  storing  the  resulting  tags  in  a  separate 
file.  New  columns  of  data  are  retrieved  until  the  entire 
large  matrix  of  data  has  been  processed. 

Figure  4  is  an  example  of  a  set  of  elevation  values.  Figure 

5  is  its  corresponding  matrix  of  tagging  values. 
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Figure  4 
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Figure  5 


The  algorithm  requires  th*t  each  of  the  values,  to  be  compared, 
has  a  full  eight  point  surrounding  neighborhood.  Therefore, 
the  first  and  last  column  and  row  of  data  are  not  processed. 
For  example,  the  six  row  by  seven  column  elevation  matrix  of 
Figure  4  has  a  corresponding  tagging  matrix  comprised  of  four 
rows  and  five  columns  (Figure  5). 

When  the  resulting  matrix  of  tagging  values  is  examined,  both 
the  location  of  the  gullies  and  their  directional  orientation 
can  be  plotted.  It  must  be  noted  that  the  directional  flow 
of  the  stream  channels  is  approximately  90  degrees  from  the 
position  of  the  neighborhood  which  triggers  the  tag.  This 
idea  is  demonstrated  in  the  above  figures.  Consider  north 
to  be  at  the  top  of  the  page.  From  the  tagging  rules,  the 
value  0010  indicates  that  the  neighborhood  north  and  south 
of  the  value  triggered  the  tag  therefore,  the  channel  would 
have  an  east-west  orientation.  If  similiar  tagging  values 
are  connected,  then  the  orientation  trend  of  the  gullies  can 
be  plotted.  With  a  few  modifications  to  the  algorithm,  a 
similiar  ridgeline  orientation  map  could  be  produced. 


TEST  DATA 

In  the  testing  of  this  algorithm,  three  elevation  data  sets 
were  used.  Two  were  manually  created  by  overlaying  a  map 
with  a  transparent  grid  and  interpolating  the  elevations 
directly  from  the  contour  lines  at  the  grid  intersection 
points.  The  grid  had  a  fineness  of  100  points  per  square 
inch.  The  third  data  set  was  a  subset  of  the  United  States 
Geological  Survey's  Digital  Elevation  Model  for  the  Koopa, 
California  area. 

The  first  test  case  was  a  matrix  of  900  elevation  values  (30 
by  30),  which  was  developed  from  a  map  of  the  Fort  Belvoir, 
Virginia  area.  This  map,  Fairfax  County/Section  108-2,  which 
was  compiled  by  Survey  and  Design,  Inc.  (Fairfax,  VA )  for 
the  Mapping  and  Graphics  Division  of  Fairfax  County  Government, 
had  a  scale  of  1  inch  to  500  feet.  Therefore,  the  elevations 
were  read  every  50  feet  for  an  area  of  1500  feet  by  1500 
feet.  The  contour  interval  was  5  feet.  This  region  was 
chosen  as  a  test  site  due  to  the  frequency  and  orientation 
of  the  drainage  channels.  One  intermittent  stream  was 
symbolised  on  the  map,  but  other  gullies  can  be  inferred  from 
the  V  indentations  of  the  contour  lines  (Strahler,  1975). 
This  particular  area  was  a  good  test  site  because  the 
delineation  of  various  channel  orientations  was  tested. 

The  second  test  area  was  a  1200  value  elevation  matrix  (30 
by  40)  which  was  devaloped  from  a  1:50,000  scale  Defense 
Mapping  Agency  Topographic  Center  sheet  of  None,  Alaska 
(Series  Q701  571  Sheet  1650  II  -  Cross  Country  Movement  Sheet 
of  the  Experimental  Alaska  Terrain  Study).  The  same  transparent 
10  line  per  inch  grid  was  used.  Therefore,  the  1200  values 
were  read  from  a  12,500  by  16,667  foot  area  with  data  points 
read  every  417  feet.  The  contour  interval  was  50  feet.  This 
area  was  chosen  because  the  stream  ehannela  tended  to  be 
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^engthy  and  have  a  variety  of  flow  orientations. 

The  third  test  case  was  a  12,000  element  matrix  (150  by  80) 
which  was  subsetted  from  the  large  USC5  Digital  Elevation 
Model  for  the  area  of  tioopa,  California.  The  spacing  between 
the  elevation  values  was  300  feet.  This  test  case  was  used 
to  demonstrate  that  automatically  produced  elevation  data 
sets  could  be  handled. 


DISCUSSIOII 


The  following  four  figures  exhibit  a  sample  of  the  algorithm 
results  using  data  from  the  second  test  case.  Figure  6  is 
a  small  area  of  the  Nome,  Alaska  1:50,000  scale  map.  A 
section  of  the  test  region  is  outlined.  The  corresponding 
elevation  data  and  resulting  tag  matrix  are  Figures  7  and  8 
respectively.  In  Figure  9,  the  locations  and  orientations 
of  the  various  gullies  are  plotted  on  the  map  from  the 
information  given  by  the  tag  matrix.  The  algorithm  automatically 
tagged  both  the  gullies  which  were  symbolised  on  the  map  and 
those  that  could  be  inferred  by  the  contour  lines.  Similiar 
results  were  demonstrated  in  the  other  test  cases. 
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Figure  7 
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Figure  8 


Figure  9 


Currently,  the  algorithm  is  limited  to  those  gullies  not 
exceeding  the  spacing  of  the  elevation  values.  Wide  streams 
channels  will  be  omitted.  For  example,  if  the  elevation 
values  are  taken  at  a  spacing  shown  by  the  x's  in  Figure  10, 
and  the  stream  has  the  following  profile,  then  it  will  be 
missed  by  the  program.  Further  algorithm  sophistication  is 
necessary  to  correct  this  problem. 
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Figure  10 


Figures  6  through  9  demonstrate  another  limitation  of  the 
current  algorithm  From  the  results,  it  appears  that  in  the 
western  section  of  the  test  area,  there  is  a  stream  which  is 
oriented  north-west  to  south-east.  Actually,  there  are  two 
streams  which  are  separated  by  a  ridgeline.  Therefore,  it 
is  planned  to  incorporate  the  ridgeline  data  with  the  gully 
location  and  orientation  tags  ~before  plotting  the  results. 
Uith  the  ridgeline  information  added,  watershed  zones  could 
also  be  delineated. 

As  with  most  work  with  elevation  matrices,  two  of  the  important 
factors  in  this  channe 1 / r idge  1  ine  delineation  technique  are 
the  quality  of  the  elevation  data  and  the  distance  separating 
the  individual  elevation  values.  If  the  input  elevation  data 
has  errors,  then  the  resulting  gully  map  will  not  correspond 
to  the  actual  ground  it  represents.  Similiarly,  if  the  area 
is  un de r s amp  1  e d  ,  then  the  resulting  map  may  miss  actual 
gullies  and/or  ridges;  since,  the  increased  distance  between 
elevation  values  results  in  a  generalization  of  the  topography. 


CONCLUSIONS 

From  these  preliminary  results  of  the  previously  discussed 
test  sets,  it  can  be  concluded  that  the  process  of  automatically 
delineating  gullies  solely  from  elevation  data  has  been 
successfully  begun.  Further  sophistication  of  this  gully 
delineation  algorithm  is  planned,  including  the  incorporation 
of  ridgeline  data  and  wide  channel  searches. 

In  the  future,  the  automated  delineation  of  a  region's  drainage 
channels  and  ridgelines  will  reduce  labor  intensive  efforts, 
vhich  occur  in  today's  mapping  operations. 
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